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^ ■ Abstract 

■ We present a study of the existence, stability and bifurcation structure of families of dark 

, breathers in a one-dimensional uniform chain of spherical beads under static load. A defocus- 

o- 

,—1 ■ ing nonlinear Schrodinger equation (NLS) is derived for frequencies that are close to the edge of 

^-H ■ the phonon band and is used to construct targeted initial conditions for numerical computations. 

. ^ . Salient features of the system include the existence of large amplitude solutions that bifurcate with 

X: 

?H . the small amplitude solutions described by the NLS equation, and the presence of a nonlinear 

^ I 

L - _ ■ 

instability that, to the best of the authors knowledge, has not been observed in classical Fermi- 
Pasta-Ulam lattices. Finally, it is also demonstrated that these dark breathers can be detected in 
a physically realistic way by merely actuating the ends of an initially at rest chain of beads and 
inducing destructive interference between their signals. 
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I. INTRODUCTION 



Granular crystals, which consist of closely packed arrays of elastically interacting parti- 
cles, have sparked recent theoretical and experimental interest. On one hand, the relevant 
mathematical model is difficult to analyze due to, among other things, the lack of smooth- 
ness of the corresponding potential function, calling for new ideas to study this nonlinear 
lattice problem On the other hand, the ability to conduct experiments for a wide range 

of setups, including chains consisting of beads of various geometries and masses, nonlinear- 



ity strength, and even dimension [3 



-|6|, make the subject of granular crystals excitin g fo r 



its potential relevance in applications inc 



actuating devices 



uding shock and ener gy a bsorbing layers 



0]io 



15 



io|, 

acoustic lenses [l2|, acoustic diodes [l3|, and sound scramblers 



16| . to name a few. 



One fundamental structure known to exist in granular crystals is the so-called intrinsic 
localized mode, or discrete breather. Discrete breathers are time-periodic solutions of the 
underlying equations of motion which are localized in space. Generally speaking, they are 
well studied, and are known to exist in a host of nonlinear lattice models jx^j. The most 
commonly studied breather is one with tails decaying to zero, which is often referred to as 
a bright breather. The term "bright" is used, as the solution can be viewed as a discrete 
carrier wave with amplitudes modulated by a bright soliton. It is then natural to consider 
a discrete wave modulated by a dark soliton, which is then in turn called a dark breather, 
see Fig. [1] for an example. Although bright discrete breathers are known to exist in granular 
crystals (in dimer or higher periodic configurations, and in monomer chains with defects for 



example 



18N21|). the properties of dark discrete breathers remains an open question, and 



is the subject of this work. 

It should be highlighted here that in monomer i.e., homogeneous chains (for which, 
there is purely an acoustic band), the conditions established for the bifurcation of discrete 
breathers 17|] do not allow the bifurcation of bright breathers. For that reason the latter 
have only been identified in heterogeneous configurations with different periodicities 19|-|21|. 
or in monomer settings bearing defects 181], where the localization is not intrinsic to the 
chain but rather is supported as an impurity mode by the defect. The only possibility that 
exists in the homogeneous chain, under static load (i.e., in the linearizable limit) for intrinsic 
localization may arise in the form of the dark breathers proposed herein. For this reason, we 
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FIG. 1. Top: A site-centered dark breather solution with Ub = 1-9 in the strain variable ?/„ = 
Un+i — Un (left), renormalized strain variable /U — |y„|, where ^ is the breather amplitude (middle) 
and in displacement variables u„ (right). The markers indicate where the solution takes values on 
the lattice. The solid line is shown for clarity. Bottom: Same as top row but for the bond-centered 
solution. 



believe that such states are fundamental ones in the study of granular systems and merit an 
investigation of their existence properties, linear and nonlinear stability, as well as of their 
dynamics. These topics form the focus of the present work which is structured as follows. 
In section 2, we present the general properties of the model. In section 3, we analyze the 
model near the band edge of its linear acoustic spectrum, and derive a defocusing nonlinear 
Schrodinger equation characterizing its envelope dynamics. In section 4, we specify the finite 
dimensional lattice dynamical system to be studied and its bifurcation analysis is carried out 
in section 5. While sections 6 and 8 consider variants of the center position of the breather 
and of the associated boundary conditions, section 7 focuses on the linear stability prop- 
erties of the states. Section 9 analyzes the validity of the analytical approximation, while 
section 10 makes connections with current experimental settings, presenting a scheme for 
the potential realization of the dark breathers. Finally, section 11 presents our conclusions 
and some challenges for future work. 
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II. MODEL 



The model describing the dynamics of a one dimensional (ID) homogeneous chain of 
spherical beads is given by js], [gl : 

with, 

V^^{x) = A[6o-x]f, (2) 

where n E I, with / a countable index set, u„ = G M is the displacement of the 

n-th bead from equilibrium position at time t, A is a. material parameter (depending on the 
elastic properties of the material and the geometric characteristics of the beads |3|), M is 
the bead mass and Sq is an equilibrium displacement induced by a static load Fq = Ad^^"^. 
The bracket is defined by = max(0,a;). Equation ([T]) with / = Z has the Hamiltonian, 



H = Y^ \mvI + \/Gc(^n+l - Mn), 



2 

neN 

where Vn = Un- Considering small amplitude solutions, i.e., 

I'^n— 1 'Un\ 
^0 

implies that we can Taylor expand Vqq{x) . Keeping terms up to fourth order yields an 
approximate model, 

MUn = Vpp^{Un+l - Un) " V^p^{Un - ^n-l), (3) 



where 



and 



^Fpu(a;) = K2X + Ksx^ + K^x 



K. = \Ab]!\ K, = -Ay-'l\ K, = -A^^^-'l\ (4) 
Equation ^ is an example of the well studied Fermi-Pasta-Ulam (FPU) model. It is has 



been shown by James 



22 



23l ] (see also the discussion in [l7|) that small amplitude bright 
breathers exist in the FPU model for frequencies above the phonon band if and only if 
B := 3K2K4 — 4i^| > 0. For coefficients such that B < 0, small amplitude dark breathers 
were shown to exist for frequencies within the phonon band. The coefficients (jl]) associated 



to the granular crystal model correspond to S < 0, and thus the existence results for dark 
breathers can be carried over for sufficiently small amplitudes (see e.g. Theorem 5 of 23|). 

In this paper, we begin with small amplitude solutions, where the phenomenology of the 
FPU lattice and granular crystal lattice are very similar. However, we will see that as we 
depart from the familiar small amplitude limit, a host of new phenomena will emerge that 
are particular to the monomer granular crystal model and its rather special nonlinearity 
bearing the unusual 3/2 power, but also the tensionless characteristic (namely, that there is 
only force between the beads when they are in contact). 



III. DERIVATION OF THE DEFOCUSING NLS EQUATION IN THE STRAIN 
VARIABLE FORMULATION 

It is more natural to derive the nonlinear Schrodinger equation in terms of the strain 
variable y„ = Un+i — Un- We note that for a finite chain, a displacement variable solution 
can be recovered using the relation, 

n-l 

Un = Ui + '^yj, (5) 

i=i 

where an arbitrary choice for the first node must be made. Equation ([1]) has the following 
form in the strain variables, 

Mjjn = -V^ciVn+l) + ^V^ciVn) " ^Gc(^n-l). (6) 

The FPU equation can again be derived, yielding, 

Mjjn = ^FPu(l/n+l) - 2l^FPu(l/n) + l^j^pu (Z/n-l) • (7) 

The linear problem (i.e. when = K4 = 0) is solved by, 

for all /c G M where u and k are related through the dispersion relation. 



u{kY = 4K2/Msm\k/2), 

such that the cutoff point of the phonon band is 2y^K^jM. 
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Consider the standard NLS multiple scale ansatz, which up to first order has the form, 

Vnit) ^ Mt) ■= eMX, T)e*(^o"+"'o*) + c.c. , X = e{n + ct), T = eH, (8) 

where e ^ 1 is a small parameter, effectively parameterizing the solution amplitude (and 
also its inverse width). The substitution of this ansatz into ([7]) and equation of the various 
orders of e leads to the dispersion relation uq = u^ko), the group velocity relation c = u'^k^), 
and the nonlinear Schrodinger equation, 

tdTA{X,T) + U2d^^A{X,T) + usA{X,T)\A{X,T)\^ = 0, (9) 

where 1^2 = —oj"{ko)/2 > and 



while 



^ _ uj{ko) / uj{2ko) cu(2fco) 

~ 2 \2uj{ko) -u{2ko) 2u{ko) + uj{2ko) ^ ^ 

2u'{0) 2uj'{0) 



uj'{ko) -Lo'{0) uj'{ko) +io'{0). 
Full details of the derivation of the NLS equation starting from the homogeneous FPU 



equation, including the higher order terms of the ansatz, can be found e.g. in 



24| . Since we 



seek standing wave solutions, we choose the wavenumber to be at the edge of the phonon 



band ko = it, such that the group velocity vanishes, ujq = 2^y K2/M, and 

^3|,„=. = 3K2i^4-4K| = S. 

We already noted that B < (and hence z/3 < 0) for coefficients defined in (jl]). Thus, the 
NLS equation iQ is defocusing. Steady states of (Q have the form. 



A{X,T) = A{X)e 
where k G M and A{X) satisfies the Duffing equation. 



aii = -l(|K|i-|i.3|i'), (12) 
where k < 0. Equation ( !T2|) possesses heteroclinic solutions which have the form. 



A{X) = VK/z/gtanh ( a/ - 
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Reconstructing the original solution from our ansatz of Eq. (|8]), we have the following ap- 
proximation, 

yn{t) ^ eie'"^e'('=°"+"«*) + c.c. 



= 2£:(— 1)". / — tanh ( i / — '-e(n — Xq) ) cosfw^t) (for ko = it) 

V 1^3 V V 2z/2 / 

where Ub = Uo + Ke"^ is the frequency of the breather, k < is a fixed but arbitrary parameter 

and Xq = exq is an arbitrary spatial translation. When k < 0, the frequency of the breather 

lies within the phonon band. We note that the same approximation (up to first order) can 



be obtained by performing a center manifold reduction (see e.g. 22|, |26|] with n = —U2 



e = y/—fi and A = M = 1 such that U2 = 1/4:). It is relevant to mention that the defocusing 



NLS equation can also be derived in the displacement variable formulation [27 1. 

If z/3 > (which is not possible here), then one would need k > 0, which, in turn, would 
result in the existence of homoclinic solutions and thus of bright breathers with frequencies 
above the phonon band. 

IV. THE FINITE DIMENSIONAL SYSTEM AND ITS CONSERVATION LAWS 

The boundary conditions (BCs) of the finite dimensional system used for the numerical 
simulations will play an important role since the dark breathers have non-decaying tails. 
In terms of the strain variable, a logical choice is a periodic boundary condition yo = i/n , 
Vn+i = Vi- To check our numerical results, we also carry out simulations in the displacement 
variables. The following system of + 1 equations is consistent with the system in the strain 
variable representation and periodic boundary conditions. 

Mill = M^o - {un+1 - un))+'^ - M^o - {u2 - ^^l))+^ 
Mii, = A{6o - (u„ - u.,.i)ff - A{6o - (n„+i - Un)f/\ n e {2, N} 
Mun+1 = A{6o - {un+i - un))+^ - A{6o - {u2 - ^^l))+^ 

The redundancy of the last equation implies un+i = ui-, such that the above may be written 
as a system of N equations, 

^3/2 ./r / nn3/2 



Mill = A{6o - {ui + cit + C2- un))+'^ - A{do - {u2 - Mi))+ , 



Miln = A{6o - {Un - Un-l)n - ^(^0 " (^n+l " Mn))7 , U G {2, N - 1} 

MiiN = A{6o - {un - un-i))+'^ - A{6o - {ui + Cit + Ca - ^^^^))+^ 
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FIG. 2. Left: Renormalized l'^ norm of the site-centered solution versus ujt with a vertical center 
of c = 0. The vertical line shows the edge of the phonon band lo = 2. The upper branch of the 
curve reaches past the edge of the phonon band and increases indefinitely. Right: Same as left, 
but for the bond-centered solution. Regions where a real instability is present are indicated by a 
dashed red line (see Sec. IVIII for details). 

where ci = nAr+i(0) — ni(0) and C2 = UAr+i(0) — ui{0). This system maintains the con- 
served quantities of the infinite system. For example, if ci = 0, then the Hamiltonian and 
mechanical momentum, 



H 



N-l 

E 

n=l 



Ur. 



N 



Vgc(mi + C2 - Mat), P = y^Mf„, 



n=l 



are conserved respectively. The conservation of the Hamiltonian is connected to the temporal 
translation invariance t ^ t + 5t and the conservation of momentum is a result of the vertical 
translation invariance u = u + 6u. An interesting feature is the invariance. 



Un^<yun + cn, t^a ^^H, a = {6Q — c)/6o, 



(13) 



which holds in the infinite and finite lattice using the above mentioned boundary conditions. 
This last invariance is connected with the existence of a marginal mode and it has important 
consequences in the linear stability analysis described in Sec. IVIII 



V. BIFURCATIONS 



We compute numerically exact dark breather solutions using a Newton-type solver for 
time-periodic solutions (see e.g. [l7| for details). We carry out the computations in the strain 
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variables with periodic boundary conditions with a lattice size of = 51. The associated 
Jacobian matrix is rank deficient due to the invariances of the system, and therefore an 
additional constraint is needed. It is sufficient to fix the vertical center of the solution, 

suPfg[o.r,] yijt) + inff6[o.r,] yijt) 

where Th = 27r/uJb is the breather period. For all reported results we used the parameter 
values M = 60 = I and A = 2/3 such that uq = 2. After a rescaling of time and amplitude, 
solutions for arbitrary parameter values can be recovered. The norm and energy of the 
dark breathers will depend on the nonzero background of the solutions. Thus, to measure 



the solutions we consider the renormalization 



28 



yn = i^-\yn- c\ hnWp = ^^J^^ -\yn- c|^, 

ne/ 

where is defined as the amplitude of the breather, 

_ suptgjo^r^] yi{t) - inftg[o,T,] yi{t) 

The renormalized profiles can be seen in the middle panels of Fig. [H We use the analytical 
approximation (jH]) to seed the numerical solver. Ansatz (E]) with xq = represents a so- 
called site-centered solution (see Fig. [T] top) whereas the bond-centered solution corresponds 
to Xo = 0.5 (see Fig. [1] bottom). 

In the finite lattice, eigenvectors of the equations of motion linearized about the trivial 
solution also provide an option for an initial seed. The largest computed eigenvalue u of 
the linear spectrum serves as the numerical cutoff value of the frequency (which is slightly 
smaller than the exact value of uq but approaches it as A^ — )■ 00). The eigenvector of the 
A^ — 1 eigenvalue has the profile of a dark breather. The drawback of using the linearized 
solution is that there is no information regarding the correct amplitude. For this reason, 
we use ([8]) for the initial seed, which approximates both the shape and amplitude of the 
solution. 

In the left panel of Fig. [21 we show the renormalized P norm of the numerically exact 
site-centered dark breathers. The solutions do not exist for arbitrarily small ujb but rather, 
there is a fold at Ub ~ 1.87. However, we note here that this is a fold in the particular 
parameters which does not constitute a point of change in the spectral stability of the 
observed states. The top branch persists past the edge of the phonon band uq = 2, with 

9 
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FIG. 3. Left: A site-centered dark breather solution with cob = 1.9 and a nonzero center c = 0.13 
in the strain variable y„. Middle: The same solution but in the displacement variables Right: 
Renormalized P norm of the site-centered solution versus Wf, with a vertical center of c = 0.13. 
The vertical line shows the edge of the phonon band oj = 2. In contrast to the c = case shown in 
Fig. [21 the bottom branch terminates at a value below the edge of the phonon band. 

indefinitely increasing amplitude. For this reason, the solutions making up the top branch 
in the figure can be thought of as strongly nonlinear solutions, and those making up the 
bottom branch as weakly nonlinear solutions. A red dashed line indicates that the solution 
possesses a real instability, a blue solid line indicates the solution does not (see Sec. IVIII for 
more details). 



VI. DARK BREATHERS WITH A NONZERO CENTER 

The equations of motion for both the infinite system and the finite system with periodic 
BCs are invariant under the transformation 

yn-^Vna + c, t-^a'~^^H, a = {5Q-c)/5o, (14) 



where (IMj) is merely the strain variable equivalent of f[T3|) . Thus, there is an entire family 
of solutions associated with the vertical shift c. We note that a similar transformation 

for the case of FPU lattices. Figure [3] shows an example of a 



has been introduced in 



dark breather with c = 0.13 in the strain (left) and displacement variables (middle). The 
bifurcation diagram (right) is similar to its c = counterpart, but is shifted in the uj axis 
where the cutoff point is ojqq^^^ = 2a^^^. The invariance ( fT^ can also be used as another 
accuracy check for the numerical simulations. We can choose the value of the constraint 
c in the Newton method and test if it converges to the appropriate solution as predicted 
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FIG. 4. Floquet multipliers of the strongly (left) and weakly (right) nonlinear solutions with c = 
versus frequency oj for the site-centered (top) and bond-centered (bottom) types. The Floquet 
multipliers for Uh = 1.95 in the complex plane are shown in the respective insets. 



by (fT4l) . assuming we have computed the c = solution. Naturally, the stability of the two 
solutions also has the same structural characteristics, as regards intervals with real Floquet 
multipliers. 

Note that a constant vertical shift y — )■ y + c is the same as leaving y unchanged and 
changing the precompression 60^60 — c. Thus, all forthcoming discussions of changes in 
vertical center can be thought of as changes in precompression. For this reason it is clear 
that both the strongly and weakly nonlinear branch go toward the zero solution as the 
precompression goes to zero, since a — )■ as c — )■ (5o. 
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VII. LINEAR STABILITY ANALYSIS 



There is a rich theory for the stabihty properties of breathers in Hamiltonian systems 



17 



30 



3l|. Here, hnear stabihty is determined in the standard way: A perturbation Vn{t) 
is added to a breather solution and ([6]) is used to derive an equation describing the evolution 
of V. Keeping only linear terms in V will result in a Hills' equation with a time periodic 
coefficient of period Tb. Thus, the eigenvalues (or Floquet multipliers) of the associated 
variational matrix at t = determine the linear stability of the breather solution. Due to 
the Hamiltonian structure of the system, all Floquet multipliers must lie on the unit circle 
for the solution to be (marginally) stable, otherwise, the solution is unstable. 

There are continuous arcs of spectrum on the unit circle (in the infinite lattice limit), 
which can be computed from the phonon band of ([6]). Although, in contrast to the bright 
breather case, these arcs are not simply exp{iujTb), where u G [— cjo,ciJo]- This is due to the 
nonzero background of the solutions we are linearizing about. The isolated (point spectrum) 
multipliers must be computed numerically. If an isolated one is real and lies off the unit circle, 
then the solution possess a so-called real- instability (see Fig. HI top left for an example). 
Oscillatory instabilities are also possible, where the multiplier has both real and imaginary 
parts (and the corresponding multipliers come in quartets). In FPU lattices, it is known 
that numerically computed eigenvalues rnay lie off the unit circle due to the truncation of 
the infinite lattice to a finite lattice [l^, Sl To determine if the instability is a finite size 
effect, a band analysis can be carried out |l7|]. Such a computation is outside the scope of 
this work, but we did carry out numerical simulations for larger values of the lattice size A^. 
We found that the oscillatory multipliers do approach the unit circle, but level off before 
reaching it. In other words, the oscillatory instabilities seem to be genuine, albeit weak. 

Examples of the numerically computed Floquet spectrum of the site-centered and bond- 
centered solutions, for both strongly and weakly nonlinear types, are shown in the insets 
of Fig. m for Ub = 1.95. The moduli of the Floquet multipliers are also shown for various 
frequencies where real and oscillatory instabilities are present. The real instabilities are the 
larger arcs in the figure. For each solution type, there is an intricate cascade of oscillatory 
instabilities, which can be best seen for the strongly nonlinear site-centered solution. The 
magnitude of the oscillatory instabilities is small relative to that of the real instabilities. 
Therefore, in the bifurcation diagram of Fig. [2] only real instabilities are indicated, although 
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FIG. 5. Manifestation of the real-instability of the strongly nonlinear bond-centered solution for 
Uh = 2.5. The left panel shows the evolution of the n = 6 node. The right panel is a contour plot 
of the time evolution of the entire solution. The color intensity corresponds to the value of the 
strain The checkered pattern is a consequence of the time periodicity of the solution. Notice 
the density dip of the solution near n = 0, and how this structure re- arranges itself in a way 
reminiscent of a moving (dark breather) pattern after t ~ 25. 



oscillatory instabilities may be present. In Fig. [5] the manifestation of an instability of the 
strongly nonlinear bond-centered solution for Ub = 2.5 is shown, which is characteristic of 
the simulations carried out for the perturbed unstable solutions. It can be clearly seen that 
the instability essentially sets the dark breather state in motion, although it does not appear 
to correspond to a genuinely traveling breather. 

SA) described above can fail to describe instabilities that 



The linear stability analysis 
grow slower than exponentially 30|. Such "nonlinear" instabilities can emerge even if the 
linear stability analysis predicts the solution to be marginally stable. They are associated 
with degenerate eigenvalues (i.e. unit Floquet multipliers), although the presence of de- 
generate eigenvalues does not, in turn, guarantee the existence of a nonlinear instability. 
For solutions with < 6o there are, up to machine precision, six degenerate eigenvalues. 
One pair, corresponding to the phase mode Vp = {y{0),y{0)}, is associated with the energy 
conservation/ time reversal invariance of the system. The second pair associated to a pair of 
degenerate eigenvalues seems to be of the form Vt = {dy/dxo,dy/dxo} which corresponds 
to continuous spatial translation invariance. This mode, called also "translational" or "pin- 
ning" mode, is of particular interest since it suggests that (exact) traveling dark breathers 
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may exist in this model. Indeed for frequencies very close to the band edge, we found sev- 
eral intermediate solutions, which seem to be spatial translations of the same solution (i.e. 
solutions with various values of Xq in Eq. (|H])). For solutions with > Sq this invariance 
is lost, and there are only four degenerate eigenvalues in the linear stability problem. 

The last pair of modes has the form Vc = {dy /dc,dy /dc} where c is the vertical center. 
To the best of the authors knowledge, these modes have never been found/examined in 
FPU lattices. Perturbing along the direction of these modes results in algebraic growth 
in the perturbation. Thus, we conjecture that there is a nonlinear instability associated 
with this marginal mode, (which in turn is connected to the shift invariance f lT^ ). The 
heuristic argument for the presence of such an instability is as follows: if one chooses a 
small perturbation of the form e„ + c, the dynamics of the solution with a vertical center 
of c may emerge, which will have a lower frequency, due to flT^ . This is precisely what 
we observe in the numerical simulations with such perturbations. For example, the weakly 
nonlinear site-centered solution at Ub = 1.97 is linearly stable. However a perturbation with 
initial value Vn{0) = .01 max(7/(0))(r„ + c), where r„ is a random perturbation centered at 
with |r„| < 1, will grow even though the LSA predicts that it should not. Moreover, 
the projection p{t) =< 2/cxact(-,t) - 2/pcrturbcd(-, t), K > grows algebraically and is several 
orders of magnitude greater than the projections to the other linearization eigenvectors. 
This phenomenon is not observed when the above perturbation has c = for the identical 
timescale. Lastly, if the invariance ( !T4l) is absent, the nonlinear instability vanishes (see 
Sec. I Villi for a relevant example). 

We note in passing, that working on the equivalent displacement representation, see Sec. 4, 
a fourth pair of Floquet multipliers is located at (1, 0) in the complex plane. This pair arises 
from the conservation of the total mechanical momentum. 

VIII. OTHER BOUNDARY CONDITIONS 

Periodic boundary conditions are advantageous since the solution profiles closely resemble 
those of the infinite system, and invariances of the infinite system carry over to the finite 
case. However, other boundary conditions are also relevant and interesting to examine from 
an experimental point of view. The free boundary conditions Uq := Ui and Mat+i := Un 
constitute such an example. In terms of strain variables, these become fixed boundary 
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FIG. 6. Left: A site-centered dark breather solution with ujb = 1-99 in the strain variable yn with 
fixed boundary conditions yo = yN+i = with N = 101. Right: The same solution but in the 
displacement variables u„. 

conditions yo = yN+i = 0. The resulting profile resembles a multi-breather, see Fig. [61 
This solution is no longer centered at zero. Unlike the periodic BC situation, one cannot 
normalize the solution to have a zero center, since the vertical shift invariance f|T^ is broken 
in the fixed BC system. Consequently, the corresponding marginal mode discussed in Sec JVIII 
vanishes, along with the eigenvector Vc and hence the nonlinear instability is absent in this 
case. Thus, choosing a strictly positive perturbation will no longer cause the integrator to 
shift towards a solution of a different frequency, since the corresponding family of vertically 
shifted solutions no longer exists. Indeed, in the simulations a perturbation of the form e^ + c 
induces the nonlinear instability in the periodic BC system, but the exact same perturbation 
in the fixed BC system does not lead to any instability and the corresponding dark breather 
is found to be structurally robust (where appropriate), in agreement with the perturbations 
of the linear stability analysis presented above. 

IX. ACCURACY OF THE NLS APPROXIMATION 

Although the NLS approximation (IT^ was sufficiently accurate to provide an initial seed 
for our numerical solver, it is relevant to investigate the extent to which the NLS approx- 
imation represents actual solutions to (I6l). Figure [7] shows a numerically exact solution in 
the strain variable (markers), the corresponding approximation (solid line) and the enve- 
lope (dashed line), which is defined by the NLS equation. From a visual inspection, the 
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FIG. 7. Left: Analytical approximation (solid line) with an envelope function that is the NLS dark 
soliton (dashed line), and a numerically exact dark breather (markers) for lji, = 1.9 and c = 0. 
Right: The amplitude of the NLS approximation (solid line), numerical breather with fixed c = 
(dashed-dot line) and numerical breather with varying c (dashed line). 



structure is quite proximal to the exact one, although the amplitude is underestimated. We 
also directly compared the amplitudes of the analytical prediction with the numerical solu- 
tions (right panel) for various Uf,, where it can be seen that NLS approximation becomes 
irrelevant as Ub moves away from the cutoff point u = 2. It can perhaps be understood at 
an intuitive level, by inspecting Figure [71 that the NLS approximation is likely to be less 
successful for the dark breather states considered herein than the more standard case of the 
bright breathers 



171]. In the latter case, the solutions bifurcate from vanishing amplitude, 
on top of a vanishing background, while for the solutions considered herein the amplitude of 
the background is finite even close to the to the cutoff limit. Hence, the cubic nonlinearity 
approximation is expected to be less adequate even in the vicinity of that limit. Moreover, 
the fold bifurcation is not predicted by the NLS approximation and there is no approxima- 
tion for the strongly nonlinear modes. Thus, the NLS approximation is only relevant for the 
weakly nonlinear solutions near the edge of the phonon band, as expected. It is interesting 
to point out that the amplitude of the dark breather with fixed c = is not the closest one 
to that of the NLS approximation. 

The NLS approximation has been rigorously justified in the context of the FPU lattice 
m |2J], and more recently for an FPU equation with an arbitrary periodic configuration in 



25|. There, it can be shown that solutions with initial data defined by the NLS equation 
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FIG. 8. Top: contour plots of the space-time evolution of an exact weakly nonlinear, bond-centered 
dark breather with a vertical shift of c = and frequency Uj, = 1.99 (left) and an initially zero 
chain actuated out of phase at the boundaries with frequency = 1.99 (right). The color intensity 
corresponds to the strain Bottom: A zoom of the corresponding top panels between times 90 
and 130. 



will remain close to an actual solution on long, but finite time scales, where the wavenumber 
k can be arbitrarily chosen. However, those results do not apply to dark breathers, the issue 
being, as highlighted also above, the nonzero background. The results presented in this 
work suggest that the justification considered in these works may be an interesting topic 
to investigate in the context of granular crystals, but also perhaps more generally (even for 
smoother potentials) for dark breather states existing on a finite amplitude background. 
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X. EXCITING THE DARK BREATHERS AND DISSIPATIVE EFFECTS 



In an experimental setting, it is difficult to impose an initial displacement and velocity 
to an arbitrary number of beads, and thus the exact solutions described above cannot be 
prescribed as an initial condition to the lattice. Hence, we have to resort to alternative 
approaches taking into consideration the somewhat limited actuation of the finite chain 
which is currently available experimentally. In that light, we propose an approach which 
relies solely on exciting (or actuating) the ends of a finite chain, and taking advantage of 
the resulting interference effect, in order to produce a time-periodic solution with a density 
dip (i.e., a dark breather). 

To verify theoretically the presence of dark breathers when only the boundaries of a chain 
are excited, we integrate Eq. numerically with a chain of beads at rest. We integrate the 
displacement variant of the equation since this corresponds directly to the physical situation. 
The numerically exact dark breathers which are in the strain variables can be converted 
by using the relation (I5l). Since Eq. ([T]) is invariant under the transformation + s, 

with s G M, we can vertically shift the resulting displacement dark breather by an amount 
s such that the center node of the solution has zero amplitude (wmiddie = 0). The actuators 
are represented by applying out of phase boundary conditions uicft(i) = ait) cos{ujbt) + s and 
Wright = O'it) cos{u!bt + tt) — s, whcrc a{t) is a slowly monotonically increasing function that 
satisfies sup a = fi, where /i is the amplitude of the exact dark breather solution of interest 
with frequency Ub- By driving the chain in this way, the (resulting plane) waves propagating 
from the left boundary and the right boundary will cancel out at the center bead resulting 
in a solution with Umiddie = 0, as desired. Thus, the resulting symmetry in the displacement 
variables (see bottom right panel of Fig. [1]) corresponds to a bond-centered solution in the 
strain- variables (see bottom left panel of Fig. [1]). Since this is an interference experiment, 
we chose an actuating frequency that lies within the phonon band such that linear waves 
will propagate through the chain. 

The left panels of Fig. [S] show contour plots of the space-time evolution of a weakly 
nonlinear, bond-centered dark breather with a vertical shift of c = in the strains. The 
amplitude in this case is /i ~ 0.8. The checkered pattern is a result of the time periodicity. 
Notice the density dip at the center of the structure. The right panel shows the contour plot 
of the space-time evolution of the actuated chain. There is a strong resemblance to a dark 
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FIG. 9. Left: The driven chain at t ~ 110 (markers) and an exact dark breather (solid hne) for 
Uf, = 1.99. Top Right: Evolution of the middle bead for the driven chain (dashed line) and exact 
dark breather (solid line). Bottom Right: The norm of the difference between the driven chain 
and exact dark breather for time points such that 2/middlo(i) is at a maximum. 

breather for time sufficiently large (here for t > 90). 

To measure the proximity of the interference to an exact dark breather solution we com- 
pared the two directly. First, we compared the difference of the time evolution of the center 
nodes ymiddic(^) of the exact and approximate breathers (top right panel of Fig. [9]). It was 
also useful to compare the entire solution at one time slice per period of oscillation (the slice 
that corresponds to a maximum of ymiddic(i) )• We used the norm to measure difference 
of the two structures (bottom right panel of Fig. [9]). The most "dark breather like" state 
during the propagation (using these measures) is shown in the left panel of Fig. [91 Here, 
the presence of the dark breather is clear, albeit slightly disturbed by presence of the linear 
waves. Essentially here, we employ linear waves to produce the destructive interference pat- 
tern associated with the dark breathers and rely on the nonlinearity to subsequently produce 
the fundamental nonlinear state associated with such a setting. 

The timescale which is required for the dark breathers to manifest via the actuation is 
crucial, as the physical system will have dissipative effects, which may be detrimental to the 
emergence of the dark breathers. Dissipation can be modeled by modifying Eq. ([1]) as 

Miln = A[So + Un-l - M„,]+^ - A[So +Un- M„+l]+^ - ^Mn- (15) 

For small values of the dissipation coefficient (i.e. for large r), dark breathers can still emerge 
(see Fig. [10] left), however for larger values, the dissipative effects are indeed detrimental to 
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the emergence of the dark breathers (see Fig. [TO] right). Thus, it is important to consider 
physical values of the parameters. The scaling, 

Un = SoUn, t = /3i, r = /3f, /3 = 

can be used to recover the physically relevant solutions, where the variables with the tilde 
represent solutions in the idealized case of (5o = M = 1 and A = 2/3 used above. For a chain 
of 316-stainless steel beads we have, 

where E = 193-10^ {N/m^) is the Young's Modulus, R = 9.525-10"^ (m) is the bead radius, 
u = 0.3 is the Poisson Ratio, p = 8027.17 {kg/im?) is the bead density, and Fq = 20 (A^) is 
a typical precompression force which results in physical parameter values of A ^ 9.7576 • 
10^ (N/m^/^), M ^ 0.0029 (kg) and 6o ^ 1.6136 ■ 10^^ (m). The scaling parameter in this 
case becomes /3 ~ 3.9533 ■ 10"^. Thus, a breather appearing at approximately t ^ 90 would 
correspond to t = 90(3 ^ 3.6 (ms) in physical time. The tested dissipation coefficients of 
f = 100 and f = 10 would become r = 100/3 ~ 4 (ms) and r = 10(3 ^ 0.4 (ms) respectively 
in the physical scaling. Thus, if one were able to carry out experiments where the dissipation 
constant is on the order of r = 4 (ms), dark breathers should be detectable. We note that the 
experimentally relevant dissipation constant in Ref. [l^ corresponds to r = 2(ms), which 
demonstrates the feasibility of detecting dark breathers in currently available experimental 
settings of granular chains. 

XI. DISCUSSION AND FUTURE DIRECTIONS 

We have demonstrated that dark breathers are not only an interesting entity of the granu- 
lar crystal model from a theoretical point of view, but are also within the reach of currently 
available experimental technology, given the ability to actuate both ends of the chain of 
beads. On the theoretical side, a useful qualitative (although only approximate) approach 
towards the understanding of such states arises through the use of the NLS reduction. On 
the other hand, their detailed numerical study illustrates the existence of a small amplitude 
branch arising from one of the linear eigenmodes of the system, but also of a large amplitude 
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FIG. 10. Top: Contour plots of the space-time evolution of the driven chain with the same 
parameters as in Fig. [8] and with a dissipation coefficient of r = 100 (left) and r = 10 (right). 
Bottom: Comparison between the driven chain (markers) and exact breather solution (solid line) 
for r = 100 (left) and r = 10 (right). Here the times were chosen to yield the closest resemblance. 

branch which cannot be captured from our analytical considerations. The destructive inter- 
ference of the signals (at the appropriate frequency inside the pass band) induced by the two 
actuators produces a waveform proximal to the exact dark breathers, and this observation is 
found to persist even in the presence of realistic values of the dissipation within the system. 

Besides the realization of such experiments which would be of particular value, our study 
also raises several interesting open problems from the mathematical/theoretical point of 
view. Among these, we highlight a band analysis of the Floquet multipliers in order to cap- 
ture the continuous spectrum of the infinite lattice problem; also, a proof of the numerically 
observed nonlinear instability, apparently present due to a non-trivial scale invariance that 
we identified; finally, the rigorous justification of the NLS equation for dark breathers and 
especially of its regime of validity (which in our case appears to be relatively limited at least 
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in a quantitative sense). 

The readiness of the solutions to move, via random perturbations of the unstable solu- 
tions, and the apparent spatial translation invariance of the system render also plausible the 
potential existence of traveling breathers. In that light, a computational study of such states 
along the lines of earlier similar studies in the realm of the discrete NLS equation 32|-|34| 



would be of particular interest. The potential emergence of dark breathers in inhomogeneous 
chains, such as dimers 
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20| or trimers 



35[ may be another issue to computationally con- 



sider in detail, since bright breathers have been identified as relatively robust nonlinear 

19-2l|- 



states in such settings 



Such studies are currently in progress and will be reported in future publications. 
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